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Abstract 

In this paper we propose a method which provides a full description of solitary wave solutions of the Schrodinger 
equation with periodically varying dispersion. This method is based on analysis and polynomial deformation of the 
spectrum of an iterative map. Using this method we discover a new family of antisymmetric bisoliton solutions. 
In addition to the fact that these solutions are of interest for nonlinear fiber optics and the theory of nonlinear 
Schrodinger equations with periodic coefficients, they have potential applications for increasing of bit-rate in high 
speed optical fiber communications. 
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1. Introduction 

Chromatic dispersion of optical fibers is one of the major factors limiting capacity of fiber communica- 
tion systems. It results in broadening of optical pulses due to difference in velocities of different spectral 
components. The dispersion management technique [T] was proposed as a method to solve the problem of 
dispersive pulse broadening. An optical fiber link with dispersion management is composed of periodically 
arranged fiber spans with alternating signs of dispersion. The fiber spans are chosen to make the cumulative 
effect of dispersion small or even zero. Full compensation of the dispersive effects is achieved by making the 
mean value of the dispersion to be zero. In this case as a pulse propagates over one period it will be spread 
and then contracted to its original shape. 

Increasing bit-rates inevitably leads to manifestation of nonlinear properties of optical fiber. As bit rates 
increase, the temporal width of an optical pulse (bit carrier) tq decreases. On the other hand the energy of the 
pulse must remain above some critical level £ > £cr to provide minimal detection bit-error at the end of the 
line, and therefore the power of the signal {P ~ S/tq) increases as the pulse narrows P > £cv/tq. It is a well 
known, experimentally verified fact that the index of refraction of the fiber core grows linearly with the power 
of the signal, n ^ uq + aP (where no is linear index of refraction and a is a coefficient of Kerr nonlinearity) . 
With the high enough power the nonlinear response of the medium becomes important. Since dispersion 
management is a purely linear approach to compensate pulse broadening, validity of this approach becomes 
less obvious in the presence of nonlinearity. It was shown that in presence of nonlinearity the compensation 
is still possible for a special class of pulses known as the dispersion managed (DM) solitons [2]. Detailed 
information about the DM solitons and their applications can be found in |3]. 
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Soon after the concept of DM solitons was introduced it became clear that both analytic representation 
and accurate numerical description of the DM soliton shape are challenging problems. The first problem 
was addressed in many papers, where approximations of the central part of the DM soliton and asymptotic 
of oscillatory tails was considered (see [4|5|6j ). However the general solution of the this problem remains 
unsolved. The problem of an accurate numerical description of the DM soliton shape is also stimulated 
considerable effort, which was solved by in [7|8|9j using slow dynamics model of DM solitons proposed 
in [TU] . Later this approach used to address behavior of DM solitons for the systems with both dispersion 
and nonlinearity compensation [1 Ij . 

The dispersion compensation while resolving the physical problem with the pulse broadening, poses an 
obstacle for the numerical computation of such solutions. Following the standard approach the description of 
soliton consists of launching a signal and waiting for all the continuous radiation to escape. For the problem 
with dispersion management this process is extremely slow and ineffective on the tails of the soliton. Authors 
of [7] proposed solving of the averaged equation for the analysis of solutions, which for the solitary waves 
reduces from an integro-differential equation to a nonlinear integral equation, solvable by a nonlinear iterative 
procedure. This allowed for a fast way of achieving precision exceeding any necessary requirements for the 
practical applications. 

Unlike the conventional NLS soliton, the dispersion managed system supports propagation of bound pairs 
of solutions, or bisolitons. Existence of bisoliton solutions was shown in the numerical investigations of [12j 
and later confirmed experimentally in [13j . It should be noted that this experimental work considered an- 
tisymmetric bisolitons and therefore we limit our investigation to that type of bisolitons. These solutions 
are interesting for the theory of nonlinear partial differential equations, and represent a new phenomena 
in nonlinear fiber optics. They illustrate a new type of solutions of nonlinear Schrodinger equation where 
dispersion coefficient is replaced with a periodic function. The potential applications for the optical commu- 
nication systems were discussed in |13| . Bisolitons can be used to introduce a three-letter alphabet ( "empty 
slot", "DM soliton" and "bisoliton"), replacing a conventional binary alphabet. This could increase a system 
capacity without a change to the physical parameters of the system. 

The shape as well as other important characteristics of bisolitons were found in |14j using a modification 
of the numerical method proposed in [7l8l9j . The equation describing averaged dynamics of the solution 
is reduced to have a single universal parameter, which incorporates all of the physical and geometrical 
characteristics of the optical system [14]. Analysis of this equation showed the existence of two branches 
of bisolitons determined by the universal parameter, with a critical bifurcation point connecting the two 
branches. 

Investigation of DM bisolitons was based on the iteration of nonlinear map. The analysis of the newly 
discovered branch of solutions was complicated by the lack of convergence of the iterative procedure. Careful 
analysis of the map showed the existence of eigenvalues outside of the unit circle, explaining the existence 
of the growing modes of the map. In this work we propose to utilize the method of polynomial spectral 
transformations to construct a contractive map, preserving the original fixed points. The application of this 
method allowed for the discovery of the previously unknown branch of bisoliton solutions. These bisolitons 
present new solutions to nonlinear Schrodinger equation with periodic dispersion. A recent paper [16j indi- 
cates a new experimental evidence of existence of this branch of bisoliton solutions. The new found bisoliton 
solutions would allow to introduce a 4-letter alphabet, further increasing the capacity of existing lines. 

2. Basic equations 

A well established model for propagation of an electromagnetic pulse through a nonlinear, dispersive 
medium is the nonlinear Schrodinger equation 

iuz + d{z)utt +"f\u\'^u ^ (1) 

where the scalar function u{z,t) is the envelope of the signal, z is the distance along the fiber, t is retarded 
time. Without loss of generality we consider a fiber link with constant nonlinearity given by 7. Function 
d{z) is the dispersion as a function of distance along the fiber. For a dispersion managed system this is a 
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piecewise constant, periodic function which can be written in the form d{z) ^ do + di if < z < L/2, and 
d{z) = do — di if L/2 < z < L, where L is a dispersion map period, and di S> do- 

A fiber hnk with dispersion management has a characteristic time scale r^juj = (_Ldi/2)^/^. The physical 
meaning of this time is the width of a pulse such that the pulse will approximately double in width as it 
propagates over half of the period. The characteristic "residual dispersion length" Zj-d = '''dm/'^o — Ldi/2do 
is the distance when a pulse with width T^m propagating over a fiber with dispersion do will approximately 
double in width. Rescaling the variables t ~ T^m^: — ^.tt^/Pu results in the equation 

i«. + ^-uii+—Uii+{2T:f — \u\'u - 0. (2) 

L di Zftj Zni 

Zni = l/{jP) is characteristic "nonlinear length" on which the change of the phase of the pulse due to fiber 
nonlinearity will be of the order 1. 

To take advantage of dispersion management in the presence of nonlinearity, the compensation must take 
place on the distances where nonlinearity effects are small and negligible compared to the local value of 
dispersion. Therefore the dispersion management works in the linear regime, and two well separated scales 
of pulse propagation can be distinguished. On the first scale local dispersion will dominate over the weak 
effects the residual dispersion and the nonlinearity. There the pulse is rapidly broadened and brought back 
to its original width as a result of the alternating signs of local dispersion. The second scale is where the 
effects of nonlinearity and residual dispersion accumulate and become important. Thus, the spectrum of the 
signal can be considered in the form 

F[u] - q{n, z) exp (^-m^l{L/2) d^ ^^^^^ j ' (3) 

with the Fourier transform F defined as F[f{t)] = J dt exTp{iujt) f (t) . Here the exponent is the rapidly 
changing phase due to the large value of local dispersion and the amplitude q is varying slowly with z due 
to residual dispersion and the nonlinear effects. Substituting this form of u into the equation ^ and taking 
an average over a map period we obtain Gabitov-Turitsyn equation describing evolution of q [lOj 

iq^ - —n^q + —R{q, fl) = 0, 
^rd ^ ^nl 

Riq, n)^ J ''"^^2^^ (^i)g(»2)g*(»3)^(r>i + r!2 - - r!)df}id02df73, 

where A = il^ + fl^ — fl^ — , and S{-) is Dirac's delta function. 

We are interested in computation of antisymmetric bisolitons which are solitary wave solutions of this 
equation. For such a solution the amplitude profile is independent of the distance along the fiber, and can be 
written as q{^, z) = </3(0)e*^^, where A is a wave number of the solitary wave. With A = l/znj the function 
If is governed by the following equation 

if + dofl'^ip = Riip, n), where do = — = f° (5) 

Thus far we have reduced this three parameter partial differential equation to an integral equation with a 
single free parameter (an alternative way to do that is presented in Appendix A). 

3. Solving the integral equation 

In this paper we study the family of antisymmetric bisolitons at different positive values of parameter do- 
In order to solve equation ([5]) we use the following iterative procedure [TS] 
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Fig. 1. Energy of the two branches of bisolitons as a function of do. The solid line corresponds to the solutions with higher 
energy (upper branch) and the dashed Une corresponds to solutions with lower energy (lower branch) 



where Q is a Petviashvilli factor [17|18j which is introduced to avoid the convergence to the trivial solution 
(fi = Qhy making that solution unstable. The operator Podd (/(^)) = (/(^) ~ /(~^)) /2 is a projection onto 
odd functions, used here to project the iterated functions onto the space of antisymmetric functions. We 
solve the equation ([6]) numerically, discretizing on a grid of M points. The iterations are stopped when 
ipn satisfies the condition 



||(p„ +dof^Vn - R{Vn,mL, < 10-^ (7) 

A Fourier transform of a sum of two Gaussian pulses, shifted apart and taken with opposite signs is used as 
a seed for the iteration procedure with do = 0.2. The value of do is decreased incrementally, with a previous 
solution used as a seed for the next value of do, until the parameter reaches zero. The same procedure was 
performed in the increasing direction from do — 0.2. We found two branches of bisoliton solutions. Figure [1] 
shows the energy of the solutions as a function of do- The solutions bifurcate near dobf ~ 0.4256. We observed 
no convergence for the values of do > do^f. 

Direct application of the scheme in equation ^ demonstrates convergence only for the lower branch of 
solutions, and does not allow to find the upper branch of solutions. In order to find the solutions from the 
upper branch we study the spectrum of the N and modify the map. 



4. Convergence of the map N 



In order to understand the convergence of the iterative procedure described in equation ^ we analyze the 
spectrum of N^^ near a fixed point. In the neighborhood of the fixed point (pfp the map A'j^ is approximated 
by the linear operator K: 

7V(^fp + *) = ^fp+X(vI')+0(v|/2). (8) 

Here the operator K lacks the complex structure and should be viewed as a linear operator over the field 
of real numbers acting on a real vector space consisting of functions Re^* and Im\l/. Let the finite set {fJ,i}, 
{\<pi)}, and {((/iil} correspond to eigenvalues, right eigenvectors and left eigenvectors of K, the discretized 
version of K. Right away we see that there is an eigenvalue 1 corresponding to a constant phase shift of 
any fixed point. This eigenvector does not effect the convergence since the equation ([7]) is unaffected by the 
constant phase shift. Therefore for convenience we introduce sets {jli} = and = {\4'i)}\\4'ti=i)j 

— {{4'i\}\{4>ti=i\ where we have removed the above eigenvector and the corresponding eigenvalue. 
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Furthermore, assume that if i < j than \fli\ > \flj\, so that p,i has maximal magnitude. The foUowing 
iterative procedure is used to obtain the N^^ eigenfunction 

i'^^) = "^^'^^ = jw:m ^ ^ 

where {4>j\ is a dual vector of \4>j), satisfying {4>j\(pk) = Sjk (here 6jk is a Kronecker symbol). We assume here 
that initial function has nonzero projection onto N*'^ eigenfunction. At each step of this procedure the 
largest — 1 modes are projected out of and in the limit as the number of iterations n goes to infinity 
the next largest mode is the only one that survives. The contribution from the M^^ vector {M > N) dies 

out as {fXMlfl'NT- 

In order to perform the procedure described in equation ^ we use the explicit expression for the operator 

k 

k{^) = l.S^fpy Re [lpI^ - L{^))\ j Iv^fppdf) + L{^) 

andi(vE')^ / '"'^^f ^ (y>fp(»i)y.fp(»2)^*(»3)+ (10) 



here dH — AQ,idVl2Anz. The operator K is not self-adjoint, therefore, its left and right eigenvectors are 
different. We obtain the the left eigenvectors of the operator K by studying the right eigenvectors of the 
operator K'^ 



C = l-Sy'fp / Re 



lHO = 

sin(A+/2) _ s{n, + n,-n,-n) (n) 
A^/2 v^fp("i)^fp("2)^ ("3) j^^^^^^ df]+ 

sin(A_/2) ^ ^ J(»i~»2-»3 + ») ,^ 

A_/2 2<^fp("i)e(»2)<Pfp(»3) dn 

2 

3 



A± = Ql±nl~nlT n^ 



Figure [5] shows how the largest magnitude eigenvalue of the operator K depends on do • Here the dashed 
line corresponds to the lower branch of solutions and the solid line corresponds to the upper branch of 
solutions. Since for all values of do the largest eigenvalue of k at the lower branch solutions is less than 1 
the map N is contractive. Therefore as long as we are close enough to the fixed point the iterations will 
converge. The situation is different for the upper branch of solutions. As clearly shown on Figure [3] for all 
values of do the largest eigenvalue is larger than 1 and therefore the map is not contractive. 

5. Construction of a contractive map 

We construct a map that retains the same fixed points as A" and at the same time a map that is contractive 
in the neighborhood of the upper branch fixed points. Consider a polynomial P{x) = Xlj=o ^j^'' ^.nd use it 
to construct a new map P: 

U (12) 

N°{ip) = LP, = A^(A^(. . . N{ip))) j times 

In the neighborhood of a fixed point this operator has the following linearization 
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Fig. 2. A plot of largest eigenvalue of X as a function of time do at the upper branch solutions (solid) and at the lower branch 

solutions (dashed). 
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Fig. 3. Plot of \p,i \ vs do for the upper branch only. The graph exhibits vertical behavior near do s» 0.4256 and near do 0.1856. 
In the first case this behavior is due to bifurcation. 
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= E + E 1^*^ = "^fp E + E E ^^^^ 1^') = 



3=0 



We choose the values of bj so that the — 1- Due to this condition the new map retains the same fixed 
point as N. Thus after n iterations we have 



^"(v'fp + ^) = ^fp + E 



(14) 
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Fig. 4. A plot of the imaginary part of eigenvalue versus the real part. Plotted first 10 eigenvalues of operator K at 'P^'^^^- 
In subfigure (a) the shaded region defines the part of space where the rest of eigenvalues lay. In subfigure (b) the real part is 
plotted on logarithmic scale. 

We choose the coefficients of the polynomial hi in such a manner that the condition max{|P(/ii)|} < 1 is 
satisfied. According to equation (|14p each mode decays exponentially, and the smaller max{|P(/ii)|} is, the 
faster the convergence. 

Our approach is a modification of the technique considered in [52], where only first order polynomials 
were needed, in our case we found that 5th order polynomials seemed to be optimal. The coefficients of the 
polynomial could be optimized for the faster convergence. This issue is discussed in the matrix Chebyshev 
approximation problem |19|20|21] . 

Figure m shows the 10 largest eigenvalues of K at the fixed point V^^^ ^ found numerically using the 
method in equation ([5]). This is a typical picture for all upper branch fixed points. This computation shows 
that the imaginary part of these eigenvalues is either very small or zero. The real part of the first 10 
eigenvalues lies between and 1.012 and all other eigenvalues lie in the shaded area. There are two key 
features of the spectrum that should be emphasized. First, jli is extremely close to 1. As shown on figure |3] 
this feature is shared by all of the upper branch solutions. Second, as demonstrated by figure[5]the following 
inequality takes place: — 1| ^ |1 — ^2!- The decay of the first mode is proportional to P(/ii)" where n is 
the number of iterations performed, therefore, to make the convergence of this mode as fast as possible we 
choose a polynomial with a derivative at (1,1) near a value of — l/(/ii — 1). Due to the first feature /ii — 1 
is very small, which implies fast growth of the polynomial to the left of the point (1, 1). However, due to 
the second feature of the spectrum such a behavior is not problematic, since there is plenty of room for the 
polynomial to return close to zero. Thus, consider the polynomial presented in the figure [6] and defined as 



i=l ' 1=1 

and let n = 0.0, ra = 0.1, rg = 0.3, r4 = 0.46, rg = 1.01 



(15) 



This is a fifth order polynomial with |P(/ii)| < 1 and P(l) = 1. Following the earlier discussion we can see 
that this polynomial does in fact has large slope at (1, 1), and yet there is enough room for the values to 
drop back down as x approaches /i2- Also note that this is only a fifth degree polynomial, and therefore one 
iteration of P(A^) requires only five iterations of N. 
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Fig. 5. A plot of fli (solid) and 112 (dash-dot) versus do- This plot indicated a large gap between the first two eigenvalues for 
the upper branch solutions. 




Fig. 6. This is a plot of a polynomial suggested in equation l|15|l . Dashed lines indicate values of 1 and 
values eigenvalues of P(A'^) will be greater than 1, in which case the new map will not be contractive. 
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To catch the solution on the upper branch we have to be as close as possible to the solution before we start 
our iterative process. Before the computing the upper branch, the lower one were obtained. Because of the 
bifurcation structure is was natural to expect that the function 2ipj;^_^^ — ip^^ could serve as a good starting 
seed for the upper branch solution at do. Even with simple iterations (without the polynomial deformation 
of the map in order to make it contractive) it was clear that there exists another branch of solutions, as for 
the first few iteration the residual ([7]) was quickly decreasing (diverging in the long run.) 

The behavior of jli illustrated in Figure [3] suggests that there is another bifurcation point with dobf ~ 0.18 
and yet another branch of antisymmetric bisoliton solutions. There is strong evidence of that from applying 
the shooting 2ip^^^^ — ipj^^ already from the computed upper branch. 
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Fig. 7. Directly solving NLS equation with initial condition mq (solid) from equation I I16I I with lower branch solution ipj^^_Q 214 
do = 0.0125, di = 1.25, L = 1.2 and 7 = 1. The result of propagation mq over 3000 periods is represented by dashed line, 
showing three orders of magnitude agreement between the initial and final pulses. Here P ~ 0.0345 so that = 45 and 

^nl = 29- 

6. Direct simulation of bisoliton propagation 

As a demonstration that the computed sohitions of the averaged equation ^ are sohtary wave approxi- 
mations of the equation ([1]) we perform a direct simulation of (HI) with the following initial condition 

(remember that dg = dQ/{'yP){diL/2)). Here is the inverse Fourier transform. The pulse is launched 
at 2 = L/4 because as equation ([3]) indicates at this point the phase of the solution tpij^ is constant. Notice 
that parameters do, di, L and 7 are determined by the physical system, while P is a free parameter whose 
value is limited from below due to the upper bound on do — dobf- 

doh{di iL 

The inequality on the right must be maintained in order to stay in the regime where Gabitov-Turitsyn 
equation is valid. Figure [7] shows the result of propagating uq over 3000 dispersion map periods. This 
distance exceeds the distances needed for information transmission. The solid line represents the initial 
pulse Mo(i) of equation (|16p launched at z — L/A. The pulse uq is constructed using a fixed point </'do=0 4 
of the map N^^^^q ^. It is used as an initial condition for direct simulation of equation ([T]) solved using a 
split-step method. The dashed line shows an amplitude profile, plotted on a logarithmic scale, of the pulse 
after the distance 3000L, showing excellent agreement between the initial and final pulses and confirming 
that uo{t) is indeed a good approximation to a solitary wave solution of equation ([1]). 

7. Conclusion 

In this paper we have presented a method to calculate solitary wave solutions of the nonlinear Schrodinger 
equation with periodic dispersion and weak nonlinearity. In case when the effect of local dispersion is much 
stronger than residual dispersion and nonlinearity the Schrodinger equation is reduced to Gabitov-Turitsyn 
equation, an integro-differential equation. This equation is satisfied by the first order term in the asymptotic 
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expansion of a solution to the Schrodinger equation with periodic dispersion. Finding solitary wave solutions 
is then reduced to finding a solution of an integral equation by finding fixed points of a derived map. 

We investigated bound pairs of solitary waves, bisolitons, and showed that system has at least two branches 
of nontrivial solutions. On one of these branches, our map is not contractive. We introduce a method 
to polynomially modify the map, obtaining new map that retains the original fixed points, while being 
contractive. Effective use of this method resulted in a previously unknown branch of bisoliton solutions. In 
order to apply this method the spectrum of the map was studied. We carefully studied the parameter space 
where such solutions are possible, and found two places where solutions appear to bifurcate. 

The method of polynomial deformation as applied to our problem can be further improved. The choice of 
the polynomials can be adapted as the knowledge about the spectrum improves in the process of iterations. 
Also, a better precision on the eigenvalues may allow one to further decrease the order of the polynomial, thus 
expediting the iteration procedure. The approach taken by this method is not limited to the antisymmetric 
solutions of our equation. Higher order solutions can be found after the appropriate modification of the map. 
Thus we have used this method to study symmetric bisolitons, with the results to be published elsewhere. 

Authors would like to thank Ildar Gabitov and Linn Mollenauer for posing the problem, and Robert Indik 
and Sasha Korotkevich for the useful discussions during the work on this problem. This work supported in 
part by the State of Arizona grant TRIF (Proposition 301), NSFgrant EMSW21-VIGRE # 0602173, and 
by the Department of Energy at Los Alamos National Laboratory under contracts DE-AC52-06NA25396 
and the DOE Office of Science Advanced Computing Research (ASCR) program in Applied Mathematical 
Sciences. 



Appendix A 

Here we present another way of reducing the equation ([1]) to a dimensionless form. In the cj-space we 
introduce a slowly varying amplitude q{z,uj) as 



u(z,Lu) = g(z,w) exp ( — iw^ / dz'(^d(z') — (rf)^ 



Here the angular brackets denote averaging over z. The slow (over many periods of the dispersion map) 
dynamics of q is described by Gabitov- Turitsyn equation [10] 

iq^iuj) - {d)uj'^q{uj) +1 J — ^^^p^(5(t^i + UJ2 - ujs - uj)S{S)q{u;i)q{LJ2)q* (l^s) = 0, 



<5 = w?+u;|-w|-w^ S{S)^lexpi-iS f dz'(d(z') - 



/L/4 

For the dispersion map d{z) = do + di if < z < L/2, and d{z) — do — di if L/2 < z < L we get 
S{S) = siTi{SLdi/A)/{6Ldi/'i). We choose •^/Ldi/2 to be a unit of time and introduce the dimensionless 
frequency ^ = y/Ldi/2uj. Then S{A) = sin(A/2)/(A/2), where A = flj + - - fl^. 

We set the residual dispersion to be equal to 1/2 by choosing the unit of propagation distance: Z = 
{Ado/ Ldi)z. The nonlinearity coefficient 7 can be eliminated by choosing the unit of the pulse amplitude u 
(or q): q{z,io) = ^2dohQ{Z,n). 

The dimensionless form of the Gabitov- Turitsyn equation becomes 

iQzm l^fQin) + f ^^h^h^sin^ + ^2- nf^^Q{n,)Q{n2)Q*in,) = o. (i?) 

J {2Trr A/2 

The stationary shape solutions have the form Q{Z, 17) — e'^^ip{iVj — there is a family of solutions parame- 
terized by wave number A. 
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In the main point of this addendum is to present another way to bring the equation ([T]) to the dimensionless 
form. Here the final step in reducing the number of parameters is done by choosing the wave number of the 
inverse solution as the unit of propagation distance. This way the only remaining parameter is dimensionless 
residual dispersion do that is related to the wave number of the solution of (|17p A as do = 1/2A. The 
advantage of this approach is that we now study the family of solutions of one equation ^17\\ instead of 
considering a family of equations (parametrized by do), where for each of these equations a solution is found 
with a certain wave number. 
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